library(tidyverse)
library(lubridate)
library(dplyr)
library(leaflet)
library(plotly)
library(viridis)
theme_set(theme_minimal() + theme(legend.position = "bottom"))
options(
ggplot2.continuous.colour = "viridis",
ggplot2.continuous.fill = "viridis"
)
scale_colour_discrete = scale_color_viridis_d
scale_fill_discrete = scale_fill_viridis_dcity_files = list.files("30_cities_data")
onehundrd_city_files = list.files("100_cities_data")Time period we are interested in
period_18 = interval(ymd("2018-02-01"), ymd("2018-04-30"))
period_19 = interval(ymd("2019-02-01"), ymd("2019-04-30"))
period_20 = interval(ymd("2020-02-01"), ymd("2020-04-30"))
period_21 = interval(ymd("2021-02-01"), ymd("2021-04-30"))The daily mean PM2.5 AQI from Feb to Aprl of year 2019 and year 2020 in each city.
city_period_diff =
tibble(city = character(),
pm25 = numeric(),
pm10 = numeric(),
o3 = numeric(),
no2 = numeric(),
so2 = numeric(),
co = numeric(),
)for (city_file in city_files) {
#print(city_file)
path = str_c("30_cities_data/", city_file)
city = strsplit(city_file, split = '-')[[1]][1]
cityAir = read_csv(path) %>%
mutate(date = as.Date(date, "%Y/%m/%d")) %>%
arrange(date)
cityAir_19 = cityAir %>%
filter(date %within% period_19)
cityAir_20 = cityAir %>%
filter(date %within% period_20)
pm25_19 = mean(cityAir_19$pm25, na.rm = T)
pm25_20 = mean(cityAir_20$pm25, na.rm = T)
pm25d = pm25_20 - pm25_19
pm10_19 = mean(cityAir_19$pm10, na.rm = T)
pm10_20 = mean(cityAir_20$pm10, na.rm = T)
pm10d = pm10_20 - pm10_19
o3_19 = mean(cityAir_19$o3, na.rm = T)
o3_20 = mean(cityAir_20$o3, na.rm = T)
o3d = o3_20 - o3_19
no2_19 = mean(cityAir_19$no2, na.rm = T)
no2_20 = mean(cityAir_20$no2, na.rm = T)
no2d = no2_20 - no2_19
so2_19 = mean(cityAir_19$so2, na.rm = T)
so2_20 = mean(cityAir_20$so2, na.rm = T)
so2d = so2_20 - so2_19
co_19 = mean(cityAir_19$co, na.rm = T)
co_20 = mean(cityAir_20$co, na.rm = T)
cod = co_20 - co_19
city_period_diff =
city_period_diff %>%
add_row(city = city,
pm25 = pm25d,
pm10 = pm10d,
o3 = o3d,
no2 = no2d,
so2 = so2d,
co = cod)
}city_period_diff =
city_period_diff %>%
mutate(
city = paste(
toupper(substring(city, 1, 1)),
substring(city, 2),
sep = ""),
city25 = fct_reorder(city, pm25, .desc = T),
city10 = fct_reorder(city, pm10, .desc = T),
cityo3 = fct_reorder(city, o3, .desc = T),
cityno2 = fct_reorder(city, no2, .desc = T),
cityso2 = fct_reorder(city, so2, .desc = T),
cityco = fct_reorder(city, co, .desc = T)
)city_period_diff %>%
plot_ly(x = ~pm25, y = ~city25, type = "bar", color = ~city25,
colors = viridis_pal(option = "D")(3), visible = T) %>%
layout(title = "Feb-Aprl Daily mean PM2.5 AQI Difference, 2020 minus 2019",
<<<<<<< HEAD
=======
<<<<<<< HEAD
>>>>>>> 4a484cce9de994b0c706b861618dccfbe37f8108
xaxis = list(title = "y"),
yaxis = list(title = "x"),
barmode = 'overlay',
yaxis = list(title = "y"),
xaxis = list(title = "x"),
updatemenus = list(
list(
y = 0.8,
buttons = list(
list(method = "restyle",
args = list("x", list(~pm25)),
label = "pm25"),
list(method = "restyle",
args = list("x", list(~no2)),
label = "no2")))
))Now we will see how the distribution of daily AQI differ between time period 2019 Feb-Aprl and 2020 Feb-Aprl.
city_PM25_Distribution = tibble()
for (city_file in city_files) {
path = str_c("30_cities_data/", city_file)
city = strsplit(city_file, split = '-')[[1]][1]
cityAir = read_csv(path) %>%
mutate(date = as.Date(date, "%Y/%m/%d")) %>%
arrange(date)
city_19 = cityAir %>%
filter(date %within% period_19) %>%
mutate(period = "2019Feb-Aprl",
day = format(date,"%m-%d"),
city = city) %>%
relocate(city, period, day)
#add a fake date "2019-02-29" with all AQI values as NA
city_19 =
city_19 %>%
add_row(city = city,
period = "2019Feb-Aprl",
day = "02-29") %>%
mutate(day = as.factor(day))
city_20 = cityAir %>%
filter(date %within% period_20) %>%
mutate(period = "2020Feb-Aprl",
day = format(date,"%m-%d"),
day = as.factor(day),
city = city) %>%
relocate(city, period, day)
city_PM25_Distribution = rbind(city_PM25_Distribution, city_19)
city_PM25_Distribution = rbind(city_PM25_Distribution, city_20)
}city_PM25_Distribution =
city_PM25_Distribution %>%
mutate(period = factor(period, levels = c("2020Feb-Aprl", "2019Feb-Aprl")),
city = paste(
toupper(substring(city, 1, 1)),
substring(city, 2),
sep = ""))#This way is stupid but it works! I tried other methods but they just don't work as expected!!!
city_PM25_Distribution %>%
plot_ly(
y = ~city, x = ~pm25, color = ~period, type = "box",
colors = c(rgb(0.2, 0.6, 0.8, 0.6), rgb(0.8, 0.2, 0.2, 0.6))) %>%
add_trace(
x = ~pm10, visible = F) %>%
add_trace(
x = ~o3, visible = F) %>%
add_trace(
x = ~no2, visible = F) %>%
add_trace(
x = ~so2, visible = F) %>%
add_trace(
x = ~co, visible = F) %>%
layout(title = "Daily PM25 AQI Distribution, 2019 and 2020 Feb-Aprl",
xaxis = list(title = "Daily AQI"),
boxmode = "group",
updatemenus = list(
list(
y = 1.1,
buttons = list(
list(label = "PM25",
method = "update",
args = list(list(visible = c(T,T, F,F, F,F, F,F, F,F, F,F)),
list(title = "Daily PM25 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
list(label = "PM10",
method = "update",
args = list(list(visible = c(F,F, T,T, F,F, F,F, F,F, F,F)),
list(title = "Daily PM10 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
list(label = "O3",
method = "update",
args = list(list(visible = c(F,F, F,F, T,T, F,F, F,F, F,F)),
list(title = "Daily O3 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
list(label = "no2",
method = "update",
args = list(list(visible = c(F,F, F,F, F,F, T,T, F,F, F,F)),
list(title = "Daily NO2 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
list(label = "so2",
method = "update",
args = list(list(visible = c(F,F, F,F, F,F, F,F, T,T, F,F)),
list(title = "Daily SO2 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
list(label = "co",
method = "update",
args = list(list(visible = c(F,F, F,F, F,F, F,F, F,F, T,T)),
list(title = "Daily CO AQI Distribution, 2019 and 2020 Feb-Aprl")))
))
))## Warning: Ignoring 60 observations## Warning: Ignoring 88 observations## Warning: Ignoring 433 observations## Warning: Ignoring 60 observations## Warning: Ignoring 120 observations## Warning: Ignoring 330 observations## Warning: 'layout' objects don't have these attributes: 'boxmode'
## Valid attributes include:
## '_deprecated', 'activeshape', 'annotations', 'autosize', 'autotypenumbers', 'calendar', 'clickmode', 'coloraxis', 'colorscale', 'colorway', 'computed', 'datarevision', 'dragmode', 'editrevision', 'editType', 'font', 'geo', 'grid', 'height', 'hidesources', 'hoverdistance', 'hoverlabel', 'hovermode', 'images', 'legend', 'mapbox', 'margin', 'meta', 'metasrc', 'modebar', 'newshape', 'paper_bgcolor', 'plot_bgcolor', 'polar', 'scene', 'selectdirection', 'selectionrevision', 'separators', 'shapes', 'showlegend', 'sliders', 'spikedistance', 'template', 'ternary', 'title', 'transition', 'uirevision', 'uniformtext', 'updatemenus', 'width', 'xaxis', 'yaxis', 'barmode', 'bargap', 'mapType'